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ABSTRACT 

The  ability  to  accurately  and  efficiently  predict  the  occurrence  and  severity  of  dynamic  stall  remains  a  ma¬ 
jor  roadblock  in  the  design  and  analysis  of  conventional  rotors  as  well  as  new  concepts  for  future  vertical 
lift.  Several  approaches  to  reduce  the  cost  of  these  dynamic  stall  simulations  for  airfoils  and  finite  wings 
are  investigated.  Temporal  error  controllers,  variable  time  step  sizes,  and  feature-based  near-body  mesh 
adaptation  are  evaluated  for  their  ability  to  more  cost-effectively  predict  dynamic  stall  on  three  different  con¬ 
figurations.  A  fourth-order  temporal  controller  has  been  observed  to  provide  a  balanced  cost-accuracy  ratio, 
as  a  maximum  of  three  to  four  orders  of  magnitude  convergence  of  the  Newton  subiterations  is  obtained 
during  much  of  the  dynamic  stall  cycle.  Larger  times  steps  can  be  applied,  in  particular  during  the  attached 
upstroke  portion  of  the  dynamic  stall  cycle  with  fourth-order  temporal  convergence.  Mesh  reductions  via  a 
feature-based  two-level  adaptation  provided  a  50%  reduction  in  computational  costs  with  comparable  ac¬ 
curacy  to  a  fixed,  refined  mesh  size.  Additional  refinements  may  be  warranted  just  after  the  dynamic  stall 
onset  to  capture  the  complex  flow  field. 


NOTATION 

cioo  Free-stream  speed  of  sound  [m/s] 

AR  Blade  aspect  ratio,  b2/S 

c  Blade  chord  length  [m] 

d  Projected  diameter  of  the  airfoil  [m] 

k  Reduced  frequency,  ooc/lU^ 

M„o  Free-stream  Mach  number,  U^/a„ 

Q  Flow  variable 

R  Wing  span  [m] 

Re  Reynold  number,  U^c/v^ 

Uoo  free-stream  velocity  [m/s] 

a  Angle  of  attack  [rad] 

co  Blade  angular  velocity  [rad/s] 

1  INTRODUCTION 

Goals  for  future  vertical  lift  (FVL)  concepts  include  ro¬ 
tors  that  no  longer  encounter  dynamic  stall,  a  phe¬ 
nomenon  that  has  limited  the  forward  flight  speed  of 
helicopters.  The  ability  to  computationally  predict  the 
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onset  of  dynamic  stall,  as  well  as  its  severity,  is  nec¬ 
essary  to  design  these  new  rotor  blades.  Further¬ 
more,  the  ability  to  accurately  predict  dynamic  stall  on 
current  vehicles  W  is  also  required  as  their  mission 
roles  change  and  expand  with  new  technology  such 
as  more  powerful  engines. 

The  study  of  dynamic  stall  has  encompassed  sev¬ 
eral  decades.  Significant  progress  in  experimen¬ 
tal  analysis  was  made  in  the  since  the  1980s  by  a 
number  of  researchers,  including  but  not  limited  to 
McCroskey [2’31,  Carr[4],  Lorber[5’61  and  Piziali [71.  Re¬ 
cent  advances  in  Particle  Image  Velocimetry  (PIV) 
and  unsteady  pressure  sensors  provide  unsteady 
flow  field  and  performance  quantities  together  to  ad¬ 
vance  knowledge  in  dynamic  stall  physics  for  a  variety 
of  configurations  and  operating  conditions [8_11]. 

Many  numerical  studies  have  also  been  under¬ 
taken,  but  progress  was  limited  until  the  past  decade. 
The  advancement  of  high  performance  computing 
on  massively  parallel  processors  provides  the  abil¬ 
ity  to  simulate  these  processes  with  millions  of  de¬ 
grees  of  freedom  necessary  to  model  these  com¬ 
plex  phenomena.  These  large  resources  have  per¬ 
mitted  rotorcraft  researchers [1 2-1 6]  to  analyze  grids, 
spatial  convergence,  and  turbulence  models  for  two- 
dimensional  airfoils.  Initial  studies  focused  primarily 


on  grid  dependence  for  static  stall  with  several  turbu¬ 
lence  models  I12-13!,  as  an  initial  step  prior  to  dynamic 
stall.  Costes  et  al.  1141  revealed  non-physical  flow  phe¬ 
nomena  (chord-wise  oscillations  in  the  suction  peaks 
and  skin  friction  profiles  resembling  transition)  due  to 
numerical  errors  even  though  the  solutions  appeared 
to  be  converged. 

Subsequent  dynamic  stall  evaluations115’17’181 
identified  boundary  layer  reattachment  as  the  fea¬ 
ture  most  sensitive  to  spatial  and  temporal  resolu¬ 
tion.  Based  on  visual  inspection  of  stall  onset  and 
flow  reattachment,  they  recommended  360,000  time 
steps/cycle  x  subiterations  for  convergence.  Klein, 
Richter,  and  Altmikus1161  reaffirmed  the  sensitivity  of 
reattachment  to  time  step  size  and  recommended 
1 000  to  2000  time  steps  per  cycle  with  1 00  subitera¬ 
tions  for  temporal  convergence.  Liggett  and  Smith 1181 
conducted  a  detailed  temporal  analysis,  identifying 
a  cutoff  point  below  which  the  simulations  remained 
first  order  even  with  the  use  of  Newton  subiterations. 
They  also  identified  an  approach  to  estimate  the  ac¬ 
curacy  of  a  mesh  for  time  independence.  Burgess 
and  Jain1191  further  explored  and  confirmed  the  role 
of  subiterations  and  convergence  for  dynamic  stall. 

The  lack  of  fidelity  in  turbulence  and  transition 
closures  have  been  the  identified  as  primary  cul¬ 
prits  in  the  lack  of  success  in  prior  dynamic  stall 
simulations.  Unsteady  Reynolds-averaged  Navier- 
Stokes  (URANS)  turbulence  models  are  statistical 
closures.  Two-equation  models  have  been  shown 
by  numerous  researchers  to  be  superior  to  one- 
equation  models,  the  latest  of  these  by  Kaufmann  et 
al. [201.  Improved  correlations  with  experiment  are  ob¬ 
tained  when  the  URANS  models  are  replaced  with 
delayed  detached  eddy  simulation  (DDES)  or  hy¬ 
brid  URANS-Large  Eddy  Simulation  (LES)  turbulence 
closures118,211.  Transition  remains  elusive  and  de¬ 
pendent  on  the  configuration  and  operational  condi¬ 
tions,  with  some  researchers  reporting  improved  cor¬ 
relations  with  experiment1201  when  transition  is  ap¬ 
plied,  while  others  do  not 1221 . 

In  all  of  these  recent  dynamic  stall  simulation  ef¬ 
forts,  researchers  agree  that  an  accurate  simulation 
of  dynamic  stall  requires  very  fine  meshes  and  small 
time  step  sizes.  The  meshes  are  typically  much  more 
refined  than  industry-level  rotor  meshes,  up  to  an  or¬ 
der  of  magnitude.  The  dense  mesh  spacing  required 
in  two-dimensions  can  be  coarsened  when  three- 
dimensional  relief  effects  associated  with  separated 
flows  are  included.  However,  the  temporal  require¬ 
ments  translate  to  one  to  two  orders  of  magnitude 
smaller  physical  time  steps  than  are  usually  applied 
to  rotors,  as  well  as  an  order  of  magnitude  increase  in 
the  number  of  Newton  subiterations  in  dual  time  step¬ 
ping  schemes. 

These  spatial  and  temporal  requirements  are  a 
bottleneck  in  achieving  accurate  rotor  simulations  that 


include  dynamic  stall.  Significant  reduction  in  the 
computational  cost,  while  maintaining  the  accuracy  of 
the  aerodynamic  performance  predictions  are  sought. 
This  effort  investigates  several  different  approaches 
that  may  provide  increased  speed  and  memory  re¬ 
duction  through  the  application  of  temporal  and  spa¬ 
tial  adaptation  and  control. 


2  NUMERICAL  MODEL 

The  numerical  simulations  were  conducted  using 
OVERFLOW  2.2,  a  structured  solver  with  Chimera 
overset  grid  capabilities1231.  For  the  computations 
in  this  effort,  spatial  terms  were  discretized  using 
a  fourth-order  central  difference  algorithm  incorpo¬ 
rating  a  diagonalized  Beam-Warming  scalar  pentadi- 
agonal  scheme.  Second-order  temporal  integration 
was  achieved  by  applying  Newton  subiterations  to  a 
first-order  implicit  Euler  scheme.  Artificial  dissipa¬ 
tion  was  included  using  the  spectral-based  dissipation 
scheme. 

The  second-order  temporal  discretization  dis¬ 
cussed  above  has  been  implemented  in  a  large 
number  of  flow  solvers,  including  FUN3D1241  and 
OVERFLOW1251.  The  user  must  select  the  number 
of  subiterations  that  is  applied  at  each  physical  time 
step.  A  large  number  of  pseudo  time  steps  will  en¬ 
sure  proper  temporal  convergence,  but  will  also  be 
prohibitive  in  terms  of  CPU  cost.  It  is  therefore  of  the 
utmost  importance  to  develop  a  criterion  that  will  re¬ 
cover  the  full  temporal  accuracy  of  the  scheme  while 
maintaining  a  reasonable  cost.  The  temporal  er¬ 
ror  introduced  by  the  backward  difference  formulation 
(BDF)  schemes  can  be  estimated  by  examining  the 
residuals  obtained  with  two  different  levels  of  approx¬ 
imations  of  the  time  derivatives 1261 .  The  time  deriva¬ 
tive  in  the  Navier-Stokes  equations  can  be  written  as 

(1)(^)A  =  i(^+lQm+1  +  ^Q"  +  ^-lQ"_1 

+  <j)^_2Qn~2  +  ■■■) 

(^)B  =  ^(^+iQm+1+^Q"+^-iQ"_1 

+  tf_2  cr2+...) 

where  the  superscript  A  and  B  represents  different 
BDF  schemes  (see  Table  1).  In  the  actual  solver,  A 
will  correspond  to  the  main  temporal  scheme  and  B  to 
the  BDF  scheme  with  the  accuracy  immediately  be¬ 
low  (BDF1  for  BDF2,  BDF2  for  BDF2^,  etc.).  Sub¬ 
tracting  the  two  time  derivatives  in  Eq.  (2),  one  can 
estimate  the  temporal  error  of  the  solution  from  time 


level  n  to  m  + 1: 
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and  the  temporal  error  norm  is  given  by 


(3) 
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Most  solvers  including  OVERFLOW  apply  the  L «, 
norm  to  assess  convergence.  The  subiterations  are 
discontinued,  and  the  next  time  step  initiated  when 
the  residuals  (algebraic  error)  descend  below  a  spec¬ 
ified  fraction  of  the  temporal  error  norm  E,.  A  drop 
of  at  least  one  order  of  magnitude  is  typically  rec¬ 
ommended  in  the  literature1251,  but  there  have  been 
few  investigations  to  date  that  have  quantitatively  con¬ 
firmed  this  value,  in  particular  for  dynamic  stall. 


3  RESULTS  AND  DISCUSSION 

Three  different  configurations  were  evaluated  nu¬ 
merically  in  this  study.  Two-dimensional  computa¬ 
tions  for  the  NACA0012  and  VR7  (no  tab)  airfoils  to 
assess  and  define  the  analysis  before  most  costly 
three-dimensional  analyses  were  undertaken.  Cor¬ 
relations  were  made  with  experimental  data  from 
McCroskey121  and  Kearney191,  respectively.  Three- 
dimensional  computations  to  confirm  the  observa¬ 
tions  and  demonstrate  the  methods  were  obtained  for 
the  OA209  finite  wing [111.  As  noted  previously  in  the 
Introduction  of  this  paper,  while  the  URANS  turbu¬ 
lence  closures  and  two-dimensional  simulations  are 
not  exact,  they  can  provide  insight  into  the  behav¬ 
ior  of  the  simulations  in  three-dimensions,  as  demon¬ 
strated  previously  by  Liggett  and  Smith 1181  at  a  lower 
cost  than  with  advanced  turbulence  closures. 


3.1  Temporal  Analysis 

In  this  effort,  two  types  of  temporal  adaptation  are 
evaluated:  a  temporal  error  controller  for  Newton 
subiterations  and  change  of  the  physical  time  step 
based  on  the  flow  field  and/or  simulation  features. 
The  convergence  of  the  Newton  subiterations  was 
measured  by  taking  the  norm  of  the  right  hand  side 
residuals. 

First,  the  influence  of  controlling  residual  error  dur¬ 
ing  the  Newton  subiterations  were  investigated  for 
a  NACA0012  airfoil  undergoing  dynamic  stall  at  a 
Moo  =  0.291  and  a  chord-based  Reynolds  number  of 
Rec  =  3.76  x  106.  The  angle  of  attack  motion  was  a  = 


5°  ±  10°  *  sin{2(ot)  at  a  reduced  frequency  of  k  =  0.1. 
Two-dimensional  simulations  with  the  Menter  Shear 
Stress  Transport  (SST)  turbulence  model1271  were 
performed  for  9,000  physical  time  steps  per  cycle  with 
20  Newton  subiterations,  based  on  the  recommenda¬ 
tions  of  Liggett  and  Smith 1181.  The  NACA  0012  O-grid 
consisted  of  200  points  normal  to  the  airfoil  and  971 
circumferential  points,  with  20  of  those  being  located 
along  the  blunt  trailing  edge.  The  stream-wise  points 
were  distributed  equally  over  the  upper  and  lower  sur¬ 
faces  of  the  airfoil.  The  initial  cell  spacing  at  the  wall 
was  chosen  to  ensure  thaty+  <  1  and  that  at  least  35- 
50  normal  cells  resolve  the  boundary  layer.  The  con¬ 
vective  and  viscous  terms  were  both  discretized  us¬ 
ing  fourth-order  central  differences  with  TLNS3D  dis¬ 
sipation,  while  time  derivatives  were  discretized  using 
the  EDF2opt  scheme.  The  linear  system  of  equations 
was  solved  implicitly  using  the  approach  of  Beam  & 
Warming.  The  order  of  magnitude  was  evaluated  for 
convergence  for  1-6  orders  of  magnitude,  which  was 
the  maximum  that  could  be  reached  successfully  with 
these  simulation  conditions. 

Figure  1  illustrates  the  airfoil  performance  pre¬ 
dicted  by  these  simulations  compared  with  experi¬ 
mental  data.  As  typically  observed  for  URANS  two- 
dimensional  simulations,  the  nonlinear  lift  is  over  pre¬ 
dicted  and  dynamic  stall  occurs  two  to  four  degrees 
higher  than  experiment.  The  fourth-  and  sixth-order 
error  controllers  provide  the  most  accurate  predic¬ 
tions  for  all  three  integrated  quantities  given  the  sim¬ 
ulation  and  mesh  limitations.  In  particular,  sufficient 
temporal  resolution  is  necessary  to  minimize  the  on¬ 
set  of  nonlinear  lift  and  stall,  improving  the  simula¬ 
tion  by  almost  two  degrees  angle  of  attack.  Fourth-  or 
sixth-order  convergence  is  also  necessary  to  correctly 
capture  the  reattachment  and  linear  lift  curve  slope 
on  the  upstroke.  Drag  and  pitching  moment  upstroke 
predictions  in  the  linear  region  is  not  as  sensitive  to 
residual  magnitude.  The  magnitude  of  the  negative 
pitching  moment  is  not  influenced  significantly,  indi¬ 
cating  that  mesh  refinement  or  turbulence  modeling 
plays  a  more  key  role  for  this  variable. 

The  order  of  magnitude  convergence  that  was 
achieved  during  the  final  cycle  of  the  analysis  is  il¬ 
lustrated  in  Fig.  2.  It  is  clear  from  this  figure  why 
the  fourth-  and  sixth-order  results  are  so  similar;  over 
the  portion  of  the  simulation  which  was  not  in  the  lin¬ 
ear  upstroke,  the  full  number  of  subiterations  were 
completed  before  the  prescribed  order  of  magnitude 
convergence  was  reached.  To  determine  if  the  solu¬ 
tion  accuracy  could  be  improved  by  additional  subit¬ 
erations,  the  number  of  Newton  subiterations  was  in¬ 
creased  from  40  to  1 00  for  fourth-order  convergence. 
The  differences  observed  between  the  two  simula¬ 
tions  lie  primarily  in  the  downstroke  during  separated 
flow.  The  initial  stall  event  occurred  less  than  0.5°  ear¬ 
lier,  with  a  reduction  of  about  0.1  in  lift  coefficient  at 


Scheme 
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Table  1 :  Coefficients  for  standard  BDF  schemes  up  to  third  order 


the  peak.  Similar  behavior  is  found  for  the  drag  coeffi¬ 
cient  and  pitching  moment  coefficient.  The  additional 
cost-to-accuracy  ratio  for  the  simulation  should  be  ad¬ 
dressed  in  three-dimensional  simulations. 

It  is  clear  that  first-  and  second-order  convergence 
is  reached  throughout  the  angle  of  attack  cycle  except 
the  initial  dynamic  stall  region  in  Fig.  2.  It  is  important 
to  note  that  even  when  there  appears  to  be  coinci¬ 
dental  convergence  characteristics  (e.g.,  the  fourth- 
and  sixth-order  curves  in  Fig.  2),  the  integrated  perfor¬ 
mance  characteristics  still  include  minor  differences  in 
the  region  where  the  convergence  appears  to  be  iden¬ 
tical  (Fig  1 ).  This  implies  that  the  higher  convergence 
obtained  in  other  regions  of  the  dynamic  stall  region 
may  influence  the  temporal  integration. 

The  convergence  has  a  distinctive  impact  on  the 
flow  field  behavior,  as  illustrated  in  Fig.  3  for  the 
non  dimensional  velocity  (U /Uo°)  for  the  first  (10M) 
and  fourth  (40M)  order  of  magnitude  residual  conver¬ 
gence.  As  the  airfoil  approaches  17°,  the  40M  flow 
over  the  suction  side  influences  a  larger  area  near  the 
leading  edge,  producing  an  earlier  nonlinear  lift  rise  at 
20°.  By  22°  the  40M  flow  field  has  stalled,  while  the 
1 0M  flow  field  is  still  undergoing  some  vortex-induced 
lift.  During  the  next  three  degrees  of  angle  of  attack, 
the  40M  flow  field  indicates  stronger  and  more  de¬ 
veloped  flow  field  features  associated  with  stall.  For 
example,  in  Fig.  3i  (1 OM),  the  development  and  shed¬ 
ding  of  the  vortex  at  the  trailing  edge  is  clearly  slower 
than  at  the  same  time  in  Fig.  31  (40M),  where  it  is  lo¬ 
cated  about  1/2  chord  downstream  from  the  trailing 
edge  and  weaker. 

While  Liggett  and  Smith  I18l  indicated  that,  in  terms 
of  the  integrated  performance  quantities,  temporal  ac¬ 
curacy  could  be  measured  by  time  steps  per  cycle  x 
Newton  subiterations,  this  point  is  revisited  with  the 
temporal  error  controller.  Evaluations  of  4,500  (with 
40  subiterations)  and  18,000  (with  10  subiterations) 
time  steps  per  cycle  were  examined.  As  the  number 
of  time  steps  per  cycle  increased,  airfoil  performance 
improvements  were  observed  with  the  lower  order 
of  magnitude  (10M  and  20M),  though  they  never 
reached  comparable  results  to  the  higher  order  of 
magnitude  results  (40M  and  60M).  No  change  in  the 
predictions  of  the  fourth  (40M)  and  sixth  (60M)  pre¬ 
dictions  were  observed  with  differing  time  step  size.  In 
all  cases,  the  60M  required  all  subiterations  and  the 
40M  required  almost  all  subiterations  (Fig.  4).  Reduc¬ 


tion  in  the  number  of  maximum  subiterations  did  not 
permit  the  solution  to  converge  sufficiently,  resulting 
in  less  accurate  simulations,  similar  to  results  shown 
by  Liggett  and  Smith [181. 

These  results  were  affirmed  by  a  series  simula¬ 
tions  performed  on  the  VR7  airfoil  without  a  trail¬ 
ing  edge  tab.  This  secondary  airfoil  configuration 
was  selected  to  ensure  that  the  observations  and 
results  observed  for  the  NACA0012  airfoil  were  ex¬ 
tensible  to  other  airfoils.  Similar  results  were  ob¬ 
tained  for  the  VR7  airfoil  undergoing  dynamic  stall 
at  a  M„  =  0.2  and  a  chord-based  Reynolds  number 
of  Rec  =  0.5  x  106.  The  angle  of  attack  motion  was 
a  =  10°  ±  10°  *sin(2cot)  for  a  reduced  frequency  of 
k  =  0.13.  Data  comparisons  were  made  with  exper¬ 
imental  data  obtained  by  Kearney  and  Glezer[91  at 
Georgia  Tech.  Their  data  was  obtained  at  a  lower 
Mach  number  (M„ 0  =  0.044,  Rec  =  0.3  x  106).  To  negate 
the  need  for  low  Mach  number  preconditioning,  the 
Mach  number  was  increased  for  the  OVERFLOW 
compressible  flow  solver.  While  there  are  some  dif¬ 
ferences  in  the  data,  the  data  trends  are  consistent. 
Five  and  ten  thousand  time  steps  per  cycle  with  40 
and  20  subiterations,  respectively  were  analyzed.  Be¬ 
tween  three  to  four  orders  of  magnitude  reduction  in 
the  right  hand  side  residuals  were  obtained,  similar  to 
the  NACA0012.  Lower  order  of  magnitude  error  con¬ 
troller  selections  resulted  in  less  accurate  solutions. 
As  these  results  are  similar  to  the  NACA0012,  they 
are  not  depicted. 

The  results  of  the  simulations  for  both  the 
NACA0012  and  the  VR7  airfoils  imply  that  the  phys¬ 
ical  time  step  size  could  be  increased  over  the  lin¬ 
ear  range,  which  may  result  in  a  significant  savings 
in  computational  costs,  in  addition  to  the  savings  by 
choosing  the  correct  temporal  controller.  To  examine 
this  further,  an  adaptive  time  step  size  for  the  VR7 
case  was  investigated.  A  simple  adaptation  based 
on  the  direction  of  the  angle  of  attack  was  evaluated. 
The  upstroke  was  modeled  at  5000  time  steps  per  cy¬ 
cle,  while  the  downstroke  included  10,000  time  steps 
per  cycle.  The  number  of  subiterations  (20)  was  kept 
at  the  smaller  time  step  size  (10,000  time  steps  per 
cycle)  throughout  the  simulation.  This  is  based  on 
the  results  of  NACA001 2  and  VR7  order  of  magnitude 
study  that  indicated  that  the  number  of  subiterations 
needed  to  reach  four  orders  of  magnitude  reduction 
were  fewer  than  20  over  this  range  (Fig.  4).  The  re- 


suits  of  the  simulation  are  presented  in  Fig.  5.  The 
simulation  results  are  comparable  for  all  approaches 
(fixed  and  adaptive  times  steps)  on  the  upstroke,  as 
indicated  by  the  coincident  lines.  The  switch  to  the 
smaller  time  step  occurred  at  20°.  By  the  third  iter¬ 
ation,  no  perceptible  change  in  the  integrated  forces 
and  moments  at  the  time  step  was  noticeable.  Fur¬ 
ther  evaluation  of  the  mass  and  momentum  flux  over 
this  change  were  comparable  between  the  fixed  and 
adaptive  time  steps,  indicating  no  loss  of  conserva¬ 
tion.  On  the  downstroke,  by  the  third  cycle,  the  inte¬ 
grated  quantities  are  observed  to  be  comparable  to 
the  1 0,000  time  step  predictions.  The  cost  of  this  sim¬ 
ulation  required  approximately  49%-50%  less  than 
that  of  the  fixed  time  step  of  10,000  steps/cycle  and 
20  Newton  subiterations. 

While  this  simple  time  step  adaption  is  sufficient 
for  this  case,  a  more  sophisticated  adaptive  parame¬ 
ter  would  be  required  for  full  rotor  simulations.  This 
parameter  would  preferably  identify  the  onset  of  dy¬ 
namic  stall,  in  particular  the  nonlinear  lift  just  before 
stall.  Based  on  recent  experimental  research  at  Geor¬ 
gia  Tech  by  Kearney  and  Glezer[9],  one  such  param¬ 
eter  is  the  vorticity  flux  (w,®)  that  shows  significant 
change  just  before  and  at  the  onset  of  dynamic  stall. 
Further  development  of  a  consistent  and  efficient  time 
adaptive  scheme  is  underway. 

3.2  Spatial  Adaptation  Feasibility  Study 

A  two-dimensional  OA209  airfoil  was  evaluated  with 
the  feature-based  mesh  adaptation  features  in  OVER¬ 
FLOW  2.2 128!.  The  mesh  adaptation  feature  sensor 
was  the  undivided  difference  of  the  flow  field  con¬ 
served  variables,  Q: 


(4) 


maxi= jfj 
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Q  ref 


The  dynamic  stall  was  evaluated  at  a  Mach  num¬ 
ber  Moo  =  0.31  with  a  corresponding  Reynolds  num¬ 
ber  per  chord  of  Rec  =  1.2  x  106.  The  airfoil  motion  of 
a  =  13°  ±7° sin(cot)  at  a  reduced  frequency  of  k  =  0.05 
has  also  been  evaluated  and  previously  published  by 
ONERA  and  DLRl17f,  as  illustrated  in  Fig.  6.  The 
experimental  data  (in  gray)  and  CFD  results  using 
the  Spalart-Allmaras  turbulence  models  indicate  dis¬ 
crepancies  that  are  normally  observed  using  URANS 
turbulence  closures  in  two-dimensional  simulations: 
pitching  moment  and  lift  overshoot  and  phase  lags  at 
the  onset  of  the  dynamic  stall,  as  well  as  varying  re¬ 
covery  rates.  To  obtain  these  results,  the  simulations 
required  360,000  combined  time  steps  and  subitera¬ 
tions  on  meshes  that  ranged  from  125,000  (Tau)  to 
400,000  (elsA)  points. 


Similar  integrated  results  were  obtained  with 
OVERFLOW  2.2  for  the  Spalart-Allmaras  turbulence 
model  for  a  mesh  of  2M  points,  as  observed  in  Fig.  7. 
The  time  stepping  was  1080  steps  per  cycle  with  20 
Newton  subiterations  for  the  Eulerian  flow  with  an  ad¬ 
ditional  4  turbulent  subiterations  at  each  Newton  sub¬ 
iteration.  This  result  formed  the  basis  for  the  mesh 
adaptation  study.  The  mesh  was  then  coarsened  to 
provide  a  basis  of  analysis  for  the  mesh  adaptation. 
The  initial  size  of  the  near  body  mesh  was  set  to  0.128 
M  (417  x  103  x  3)  with  an  off-body  Cartesian  mesh 
of  approximately  0.1  M,  as  illustrated  in  Fig.  8. 

As  the  simulation  progressed,  the  airfoil  was 
adapted  with  two  levels  of  refinement  at  every  ten 
time  steps,  as  illustrated  in  Fig.  9.  The  mesh  adap¬ 
tation  with  respect  to  the  angle  of  attack  indicates 
that  the  mesh  refinement  is  required  during  the  ini¬ 
tial  vortex  shedding  that  characterizes  the  nonlinear 
lift  increase  just  prior  to  dynamic  stall.  The  maxi¬ 
mum  mesh  refinement  occurs  during  the  minimum 
pitching  moment  where  vorticity  is  shed  and  found 
in  the  near  wake  from  both  the  leading  and  trailing 
edges  (Fig.  10).  Interestingly,  the  mesh  refinement 
is  minimal  during  the  reattachment  phase,  although 
Liggett  and  Smith [18]  noted  that  this  was  a  very  sensi¬ 
tive  area  temporally.  This  implies  that  reattachment  is 
likely  temporally  driven  rather  than  mesh  driven,  given 
a  sufficient  mesh  size.  The  integrated  airfoil  perfor¬ 
mance  (Fig.  10)  indicates  that  the  adaptive  mesh  can 
predict  comparable  results  to  the  larger  fixed  mesh, 
confirming  that  this  approach  is  sufficient.  The  goal 
in  these  computations  was  not  to  optimize  the  mesh 
adaptation  size,  but  to  illustrate  that  the  mesh  adapta¬ 
tion  could  capture  the  salient  features  of  the  dynamic 
stall. 


3.3  OA209  Finite  Wing 

To  demonstrate  the  efficacy  of  the  temporal  error 
controller  and  mesh  adaptive  approaches,  a  three- 
dimensional  finite  wing  has  been  evaluated.  The  as¬ 
pect  ratio  three  (3)  wing  is  based  on  the  OA209  air¬ 
foil  and  was  evaluated  in  a  wind  tunnel  using  Parti¬ 
cle  Image  Velocimetry  (PIV),  Laser  Doppler  Velocime- 
try  (LDV),  and  unsteady  Kulite  pressure  transducers 
for  static  and  dynamic  stall.  The  free  stream  veloc¬ 
ity  was  Voc  =  55  m/s  corresponding  to  a  Moo  =  0.16 
and  a  chord-based  Reynolds  number  of  Rec  =  lx  106. 
The  tunnel  blockage  was  determined  to  be  minimal, 
with  a  measured  tunnel  turbulence  of  0.05%.  The  dy¬ 
namic  stall  tests  were  conducted  with  free  transition 
for  a  =  17°  ±5°sin(a)t)  where  co  =  3.5  Hz  or  a  reduced 
frequency,  k  =  0.1.  Further  details  and  data  from  the 
experimental  campaign  can  be  found  in  Le  Pape  et 
al.  [11]. 

The  near-body  grids  extended  to  the  farfield  and 
did  not  utilize  the  off-body  grid  capabilites  in  OVER- 


FLOW;  therefore,  only  the  near-body  grid  adaption 
option  in  OVERFLOW  was  applied.  The  non-adapted 
meshes  were  first  assessed  for  static  conditions  to 
ensure  that  they  were  sufficient  to  capture  the  three- 
dimensionality  of  the  finite  wing.  The  final  un¬ 
adapted  mesh  size  was  approximately  26.5  million 
mesh  points.  A  coarsened  mesh  of  16  million  grid 
points  was  applied  as  the  baseline  for  adaptation. 

The  mesh  adaptation  feature  sensor  was  the  undi¬ 
vided  difference  of  the  flow  field  conserved  variables 
(Eq.  4).  To  control  the  grid  growth,  a  growth  ratio  of 
1 .2  was  applied  with  a  total  of  two  adaption  levels.  Er¬ 
ror  tolerances  were  adjusted  using  Newton’s  method 
to  meet  the  grid  constraints.  In  addition,  smoothing  of 
the  error  sensor  function  was  applied  to  minimize  dis¬ 
parate  growth  in  adjacent  regions,  as  well  as  weight¬ 
ing  of  the  function  that  applied  a  decay  away  from  the 
airfoil  to  minimize  refinement  in  wake  areas  that  would 
not  impact  the  integrated  performance  quantities.  The 
mesh  was  adapted  at  every  5  steps  with  a  total  of 
18,000  time  steps/cycle  (3,600  adaptations).  Only  a 
third-order  convergence  during  the  Newton  subiter¬ 
ations  was  achieved,  unlike  the  fourth-order  conver¬ 
gence  observed  for  the  two-dimensional  simulations. 
The  standard  Splarat-Allmaras  turbulence  model  was 
applied. 

The  mesh  adaptation  over  the  final  (third)  cycle  is 
shown  in  Fig.  11.  The  upstroke  mesh  stayed  rela¬ 
tively  constant  in  size  to  the  baseline  mesh  until  the 
nonlinear  lift  region  and  the  onset  of  dynamic  stall 
was  reached.  At  that  time,  the  mesh  size  increased 
by  approximately  two  million  mesh  points  until  reat¬ 
tachment  near  the  beginning  of  the  cycle.  An  exam¬ 
ple  of  how  the  near-body  mesh  adaption  occurs  dy¬ 
namically  can  be  observed  in  Fig.  12.  The  mesh, 
and  accompanying  snapshot  of  the  near-body  shed 
wake  vorticity,  for  13.42°  on  the  upstroke  and  20.53° 
on  the  downstroke  illustrate  the  changing  mesh.  For 
the  upstroke,  the  mesh  remains  relatively  coarse,  ex¬ 
cept  for  some  refinement  near  the  root.  During  the 
downstroke,  where  there  is  significant  vorticity  shed 
from  the  wing,  there  are  multiple  regions  of  mesh  re¬ 
finement  that  extend  the  entire  wing  span,  with  sec¬ 
ondary  refinement  levels  closer  to  the  surface. 

The  lift  and  pitching  moments  are  compared  in 
Figs.  13  and  Fig.  14,  respectively  with  experimental 
data!11].  The  upstroke  portion  of  the  simulation  for 
both  the  fixed  grid  and  the  adapted  mesh  are  virtu¬ 
ally  identical  for  both  lift  and  pitching  moment.  The 
primary  differences  between  the  two  simulations  oc¬ 
cur  during  the  onset  of  dynamic  stall,  recovering  near 
the  reattachment  of  the  dynamic  stall.  The  discrep¬ 
ancies  between  the  two  simulations  is  observed  at  all 
span  stations,  but  the  longest  range  of  mismatching 
results  occurs  at  the  50%  (r/R  =  0.5)  span  station. 
The  adapted  grid  simulation  does  not  recover  the  ac¬ 
curacy  of  the  fixed  grid  simulation  until  a  =  14°.  The 


abrupt  reattachment  observed  at  the  inboard  station 
is  mitigated  by  the  influence  of  the  tip  vortex  in  sta¬ 
tions  further  outboard,  as  observed  by  the  more  gen¬ 
tle  reattachment  during  the  downstroke.  The  fixed 
and  adapted  grid  results  are  comparable  as  early 
as  a  =  16°  -  18°,  though  the  magnitude  of  the  mis¬ 
matches  at  higher  angles  of  attack  is  larger  than  ob¬ 
served  at  the  50%  span  station.  The  increased  lift  that 
occurs  in  the  experiment  at  the  onset  of  the  dynamic 
stall  at  the  outboard  stations  is  captured,  albeit  less 
intensely,  by  the  adaptive  grid  but  not  by  the  fixed,  re¬ 
fined  grid.  The  magnitude  of  the  maximum  negative 
pitching  moment  is  more  accurately  predicted  by  the 
adaptive  mesh  at  the  50%  span  station,  but  is  less 
accurate  for  the  three  outboard  stations. 

The  cost  of  the  simulation  using  the  adaptive  grid 
was  approximately  50%  of  the  cost  of  the  simula¬ 
tion  with  the  fixed,  refined  mesh.  Further  optimiza¬ 
tion  of  the  number  of  refinement  levels  and  growth 
factors  are  warranted  to  more  accurately  predict  the 
separated  region  just  after  the  onset  of  dynamic  stall. 
The  complex  wake  region  may  require  more  rapid  cell 
growth  and  additional  levels,  while  it  is  apparent  that 
upstroke  region  will  not  benefit  from  further  refine¬ 
ment.  The  vorticity  flux  or  a  similar  parameter  such 
as  Q-criterion  may  identify  regions  where  additional 
levels  of  refinement  may  be  needed  to  maintain  accu¬ 
racy. 

4  CONCLUSION 

The  ability  to  predict  dynamic  stall  accurately  is  lim¬ 
ited,  based  on  the  fixed  temporal  and  spatial  require¬ 
ments  of  most  state-of-the-art  CFD  solvers.  A  study 
has  been  performed  to  investigate  the  ability  of  tem¬ 
poral  error  controllers  and  feature-based  mesh  adap¬ 
tation  to  mitigate  these  costs.  In  addition,  an  adaptive 
time  step  approach  is  proposed  to  further  reduce  the 
costs  of  these  dynamic  stall  simulations  in  two  and 
three  dimensions.  Based  on  the  study  performed  on 
three  different  airfoil  geometries  in  two  and  three  di¬ 
mensions,  the  following  conclusions  can  be  observed: 

1.  Temporal  error  controllers  can  be  applied  to 
maintain  accuracy  during  the  attached  portion 
of  the  dynamic  stall  without  the  need  to  apply 
the  maximum  number  of  Newton  subiterations  at 
each  time  step. 

2.  During  the  separated  portion  of  the  dynamic 
stall  only  2-3  orders  of  magnitude  convergence 
were  observed,  while  during  the  attached  flow 
region,  close  to  six  orders  of  magnitude  conver¬ 
gence  can  be  obtained.  Two  orders  of  magnitude 
convergence  is  insufficient  for  accuracy  without 
many  physical  time  steps  per  cycle.  Four  orders 
of  magnitude  convergence  are  recommended  as 


it  provides  comparable  accuracy  for  the  three 
time  step  sizes  examined. 

3.  Additional  cost  reduction  (25%)  can  be  obtained 
by  increasing  the  physical  time  step  without  in¬ 
creasing  the  number  of  Newton  subiterations  dur¬ 
ing  the  upstroke  (attached  flow)  portion  of  dy¬ 
namic  stall.  A  simple  adaptive  time  step  based  on 
the  direction  of  the  pitch  has  been  demonstrated. 

4.  A  50%  reduction  in  time  was  demonstrated  with 
two-  and  three-dimensional  feature-based  spatial 
adaptation.  A  two-level  refinement  with  a  growth 
rate  of  1 .2  was  sufficient  to  capture  the  upstroke 
and  recovery/reattachment  of  the  dynamic  stall. 
Additional  refinement  is  warranted  to  more  accu¬ 
rately  capture  the  complex  wake  during  dynamic 
stall  onset  and  largely  separately  flow  field  follow¬ 
ing. 

Simultaneous  temporal  and  spatially  adapting 
meshes  have  been  demonstrated,  resulting  in  cost 
savings  of  50%  to  70%  from  traditional  fixed  meshes 
and  time  step/Newton  subiteration  schemes.  Ad¬ 
ditional  development  and  optimization  of  these  ap¬ 
proaches  for  full  rotors  and  advanced  turbulence  clo¬ 
sures  such  as  hybrid  RANS-LES  methods  are  under¬ 
way. 
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2.5 


(c)  Pitching  moment  coefficient 


(a)  Newton  Subiterations  -  9,000  steps 


Fig.  2:  Number  of  Newton  subiterations  required  to 
reach  the  specified  order  of  magnitude  drop  in  the 
residuals. 


Fig.  1:  Two-dimensional  OVERFLOW  computa¬ 
tional  analyses  using  the  Menter  SST  model  for  an 
NACA0012  airfoil  undergoing  dynamic  stall  with  tem¬ 
poral  error  controllers.  A  time  step  of  9,000  time  steps 
per  cycle  and  40  Newton  subiterations  was  applied. 


(a)  a  =  17.42°  -  IOM 


(b)  a  =  19.76°  -  IOM 


(C)  a  =  21.79°  -IOM 


(d)  a  =  17.42° -40M 


(e)  a  =  19.76° -4 OM 


(f)  a  =  21.79°  —AOM 


(g)  a  =  23.41° -IOM 


(h)  a  =  24.49° -IOM 


(i)  a  =  24.98°  -  IOM 


Fig.  3:  Impact  of  number  of  order  of  magnitude  drop  on  the  flow  field. 


(b)  Order  of  Mag.  Drop  -  40M 

Fig.  4:  Number  of  Newton  subiterations  required  to 
reach  the  specified  order  of  magnitude  drop  in  the 
residuals. 
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(b)  Pitching  moment  coefficient 


Fig.  5:  Time  adaptation  during  dynamic  stall  for  a  VR7 
(no  tab)  airfoil  undergoing  dynamic  stall  at  k  =  0.13. 


Fig.  6:  Two-dimensional  computational  analyses  us¬ 
ing  the  Spalart-Allmaras  model  for  an  OA209  air¬ 
foil  undergoing  deep  dynamic  stall.  From  Richter  et 
al. [17]. 
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Fig.  7:  Demonstration  that  OVERFLOW  without  adap¬ 
tation  obtains  similar  results  to  Richter  et  al.  I17!. 


(a)  Full  mesh 


(b)  Airfoil  mesh  with  adaptation 


Fig.  8:  OVERFLOW  airfoil  mesh  for  the  OA209  airfoil 
with  near-body  mesh  adaptation. 


a,  deg 


Fig.  9:  Mesh  adaptation  for  the  OA209  airfoil  under¬ 
going  deep  dynamic  stall. 


(a)  Lift  coefficient 
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(b)  Pitching  moment  coefficient 


Fig.  10:  Two-dimensional  OVERFLOW  computational 
analyses  using  the  Spalart-Allmaras  model  for  an 
OA209  airfoil  undergoing  deep  dynamic  stall  with  and 
without  mesh  refinement. 
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Fig.  1 1 :  Mesh  adaptation  for  the  three-dimensional 
OA209  finite  wing  undergoing  dynamic  stall. 
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(a)  13.42°  upstroke 


(b)  20.53°  downstroke 


Fig.  12:  Dynamic  mesh  adaptation  and  near-body 
shed  wake  for  the  OA209  finite  wing  during  dynamic 
stall 


(a)  r/R  =  0.5 


(a)  r/R  =  0.5 


(b)  r/R  =  0.8 
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(d)  r/R  =  0.99 
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Fig.  13:  Lift  coefficient  comparison  of  the  refined 
mesh,  adapted  mesh,  and  experiment  f11!. 


Fig.  14:  Pitching  moment  coefficient  comparison  of 
the  refined  mesh,  adapted  mesh,  and  experiment  t11!. 


